
if ("`swdLocal'\Data"=="") {
	di as error "Please run init.do first."
	error 1
}
clear all
set more off
set maxvar 10000

*Open file for regressions
use "${gsdAnalysisOutput}/bscie_panel_indices_final", clear

*Keep only endline and lagged controls
keep if wave==2 & final_sample==1
drop if replacement == 1 & gender == 1

*Define strata
egen strata=group(gender state_old)

*Generate treatment streams
gen treatmentX=0
replace treatmentX=1 if treatment_stream==3

ren edu_level_b  edu_b

*Generate IV treatment x kcb_dist 
gen Z_treat_kcb=treatment*lg_kcb_dist
label var Z_treat_kcb "Treatment x (log) distance to KCB branch"

*Generate interaction for other geographic controls
gen treat_citydist=treatment*lg_city_dist
label var treat_citydist "Treatment x (log) distance to city center"
gen treat_roaddist=treatment*lg_road_dist
label var treat_roaddist "Treatment x (log) distance to primary road"
gen treat_gradient=treatment*gradient
label var treat_gradient "Treatment x gradient"
gen treat_alt=treatment*info_loc_gps_altitude_old
label var treat_alt "Treatment x altitude"
gen treat_conflict = treatment*death_x_prox_a300
label var treat_conflict "Treatment x conflict_affected_index(300km buffer)"
gen treat_edu = treatment*edu_b
gen treat_num = treatment*num_eval_b
gen treat_lit = treatment*lit_eval_b

* Generate sum of interaction and level
gen TplusD = treatment + lg_kcb_dist 
gen TtimesDminus1 = treatment * (lg_kcb_dist - 1)
gen TplusNogrant = treatment + treatmentX
gen TtimesNOgrantminus1 = treatment * (treatmentX - 1)

*Define geographic controls
global geo_controls lg_city_dist treat_citydist treat_roaddist lg_road_dist gradient treat_gradient death_x_prox_a300 treat_conflict 

global base_controls i.edu_b i.treat_edu i.num_eval_b i.treat_num i.lit_eval_b i.treat_lit
* 
eststo: xi: ivreg2 trust_index $geo_controls TtimesDminus1 i.strata (TplusNogrant = TplusD ) , cluster(boma_old) first savefirst savefprefix(fstage_)

********************************************************************************
*** Reduced form
eststo: xi: ivreg2 consumption_index c.lg_kcb_dist##i.treatment, cluster(boma_old) 
eststo: xi: ivreg2 consumption_index c.lg_kcb_dist##i.treatment $geo_controls  i.strata, cluster(boma_old) 
eststo: xi: ivreg2 consumption_index c.lg_kcb_dist##i.treatment $geo_controls $base_controls i.strata, cluster(boma_old) 
*** Reduced form
eststo: xi: ivreg2 trust_index c.lg_kcb_dist##i.treatment, cluster(boma_old) 
eststo: xi: ivreg2 trust_index c.lg_kcb_dist##i.treatment $geo_controls  i.strata, cluster(boma_old) 
eststo: xi: ivreg2 trust_index c.lg_kcb_dist##i.treatment $geo_controls $base_controls i.strata, cluster(boma_old) 

*** Define variables that we need to estimate the join significance of lg_kcb_dist + treatmentX

gen blub = treatmentX - treatment
gen blubIV = c.lg_kcb_dist#1.treatment - treatment

***Socioeconomic Survey outcomes ******************************************************************
clear matrix

mat A1 = J(1,18,.)
mat A2 = J(1,18,.)
mat A3 = J(1,18,.)
mat A4 = J(1,18,.)
mat stars1 = J(1,18,0)
mat stars2 = J(1,18,0)
mat stars3 = J(1,18,0)
mat stars4 = J(1,18,0)

mat N = J(1,6,.)
mat D = J(1,6,.)


global n 1
local replace replace
local append
foreach var of varlist employment_index consumption_index savings_debt_index bizskills_index {
	
	* LATE geo-controls
		eststo: xi: ivreg2 `var' lg_kcb_dist treatment $geo_controls  i.strata ( treatmentX= Z_treat_kcb) , cluster(boma_old) first
		
		local b_t1 = _b[treatmentX]
		display `b_t1'
		mat A$n[1, 1] = `b_t1'
		local p_t1=(2 * ttail(e(ardf_r), abs(_b[treatmentX]/_se[treatmentX])))
		display `p_t1'
		mat A$n[1, 2] = `p_t1'
		local b_t2 = _b[treatment]
		display `b_t2'
		mat A$n[1, 4] = `b_t2'
		local p_t2=(2 * ttail(e(ardf_r), abs(_b[treatment]/_se[treatment])))
		display `p_t2'
		mat A$n[1, 5] = `p_t2'
		local obs = `e(N)'
		mat N[1,1] = `obs'
		local fstat = `e(rkf)'
		mat D[1,1] = `fstat'
			
			preserve
			parmest, norestore
			keep if parm=="treatmentX" | parm=="treatment"
			gen family="socioeconomic_controls_t1" if parm=="treatmentX"
			replace family="socioeconomic_controls_t2" if parm=="treatment"
			replace parm="`var'_t1" if parm=="treatmentX"
			replace parm="`var'_t2" if parm=="treatment"
			gen num=$n
			save "$gsdRawTemp\1qvalue_ate_con_$n", replace
			restore
			
		eststo: xi: ivreg2 `var' lg_kcb_dist treatment $geo_controls i.strata (blub = blubIV) , cluster(boma_old) 		first
		local b_d = _b[treatment]
		display `b_d'
		mat A$n[1, 7] = `b_d'
		local p_d=(2 * ttail(e(ardf_r), abs(_b[treatment]/_se[treatment])))
		display `p_d'
		mat A$n[1, 8] = `p_d'	
			preserve
			gen parm = "`var'"
			gen p = `p_d'
			gen num=$n
			keep parm p num
			duplicates drop
			save "$gsdRawTemp\1qvalue_ate_difference2_$n", replace
			restore
			
	eststo: xi: ivreg2 `var' lg_kcb_dist treatment $geo_controls $base_controls  i.strata (treatmentX= Z_treat_kcb) , cluster(boma_old) first
	
		local b_t1 = _b[treatmentX]
		display `b_t1'
		mat A$n[1, 10] = `b_t1'
		local p_t1=(2 * ttail(e(ardf_r), abs(_b[treatmentX]/_se[treatmentX])))
		display `p_t1'
		mat A$n[1, 11] = `p_t1'
		local b_t2 = _b[treatment]
		display `b_t2'
		mat A$n[1, 13] = `b_t2'
		local p_t2=(2 * ttail(e(ardf_r), abs(_b[treatment]/_se[treatment])))
		display `p_t2'
		mat A$n[1, 14] = `p_t2'
		local obs = `e(N)'
		mat N[1,4] = `obs'
		local fstat = `e(rkf)'
		mat D[1,4] = `fstat'
			
			preserve
			parmest, norestore
			keep if parm=="treatmentX" | parm=="treatment"
			gen family="socioeconomic_controls_t1" if parm=="treatmentX"
			replace family="socioeconomic_controls_t2" if parm=="treatment"
			replace parm="`var'_t1" if parm=="treatmentX"
			replace parm="`var'_t2" if parm=="treatment"
			gen num=$n
			save "$gsdRawTemp\1qvalue_ate_geocon_$n", replace
			restore
	

		eststo: xi: ivreg2 `var' lg_kcb_dist treatment $geo_controls $base_controls   i.strata (blub = blubIV) , cluster(boma_old) first
		local b_d = _b[treatment]
		display `b_d'
		mat A$n[1, 16] = `b_d'
		local p_d=(2 * ttail(e(ardf_r), abs(_b[treatment]/_se[treatment])))
		display `p_d'
		mat A$n[1, 17] = `p_d'	
			preserve
			gen parm = "`var'"
			gen p = `p_d'
			gen num=$n
			keep parm p num
			duplicates drop
			save "$gsdRawTemp\1qvalue_ate_difference3_$n", replace
			restore
	local replace
	local append append
	global n=$n+1
}

***Psychological and Behavioral Survey outcomes ******************************************************************

mat O1 = J(1,18,.)
mat O2 = J(1,18,.)
mat O3 = J(1,18,.)
mat O4 = J(1,18,.)
mat O5 = J(1,18,.)

mat starZ1 = J(1,18,0)
mat starZ2 = J(1,18,0)
mat starZ3 = J(1,18,0)
mat starZ4 = J(1,18,0)
mat starZ5 = J(1,18,0)

global n 1
local replace replace
local append
foreach var of varlist psycho_index risk_index trust_index crime_index migration_index {
	

			
	eststo: xi: ivreg2 `var' lg_kcb_dist treatment  $geo_controls i.strata ( treatmentX= Z_treat_kcb) , cluster(boma_old)  first 
		local b_t1 = _b[treatmentX]
		display `b_t1'
		mat O$n[1, 1] = `b_t1'
		local p_t1=(2 * ttail(e(ardf_r), abs(_b[treatmentX]/_se[treatmentX])))
		display `p_t1'
		mat O$n[1, 2] = `p_t1'
		local b_t2 = _b[treatment]
		display `b_t2'
		mat O$n[1, 4] = `b_t2'
		local p_t2=(2 * ttail(e(ardf_r), abs(_b[treatment]/_se[treatment])))
		display `p_t2'
		mat O$n[1, 5] = `p_t2'
		local fstat = `e(rkf)'
		mat D[1,1] = `fstat'
			preserve
			parmest, norestore
			keep if parm=="treatmentX" | parm=="treatment"
			gen family="behavioral_controls_t1" if parm=="treatmentX"
			replace family="behavioral_controls_t2" if parm=="treatment"
			replace parm="`var'_t1" if parm=="treatmentX"
			replace parm="`var'_t2" if parm=="treatment"
			gen num=$n
			save "$gsdRawTemp\1qvalue_ate_con_psy_$n", replace
			restore
			
		eststo: xi: ivreg2 `var' lg_kcb_dist treatment $geo_controls i.strata (blub = blubIV) , cluster(boma_old) first  
		local b_d = _b[treatment]
		display `b_d'
		mat O$n[1, 7] = `b_d'
		local p_d=(2 * ttail(e(ardf_r), abs(_b[treatment]/_se[treatment])))
		display `p_d'
		mat O$n[1, 8] = `p_d'	
			preserve
			gen parm = "`var'"
			gen p = `p_d'
			gen num=$n
			keep parm p num
			duplicates drop
			save "$gsdRawTemp\1qvalue_ate_difference2_psy_$n", replace
			restore
			
	eststo: xi: ivreg2 `var' lg_kcb_dist treatment  $geo_controls $base_controls   i.strata (treatmentX= Z_treat_kcb) , cluster(boma_old) first 
		local b_t1 = _b[treatmentX]
		display `b_t1'
		mat O$n[1, 10] = `b_t1'
		local p_t1=(2 * ttail(e(ardf_r), abs(_b[treatmentX]/_se[treatmentX])))
		display `p_t1'
		mat O$n[1, 11] = `p_t1'
		local b_t2 = _b[treatment]
		display `b_t2'
		mat O$n[1, 13] = `b_t2'
		local p_t2=(2 * ttail(e(ardf_r), abs(_b[treatment]/_se[treatment])))
		display `p_t2'
		mat O$n[1, 14] = `p_t2'
		local fstat = `e(rkf)'
		mat D[1,4] = `fstat'
			preserve
			parmest, norestore
			keep if parm=="treatmentX" | parm=="treatment"
			gen family="behavioral_controls_t1" if parm=="treatmentX"
			replace family="behavioral_controls_t2" if parm=="treatment"
			replace parm="`var'_t1" if parm=="treatmentX"
			replace parm="`var'_t2" if parm=="treatment"
			gen num=$n
			save "$gsdRawTemp\1qvalue_ate_geocon_psy_$n", replace
			restore			
	
		eststo: xi: ivreg2 `var' lg_kcb_dist treatment $geo_controls $base_controls   i.strata (blub = blubIV) , cluster(boma_old) first 
		local b_d = _b[treatment]
		display `b_d'
		mat O$n[1, 16] = `b_d'
		local p_d=(2 * ttail(e(ardf_r), abs(_b[treatment]/_se[treatment])))
		display `p_d'
		mat O$n[1, 17] = `p_d'	
			preserve
			gen parm = "`var'"
			gen p = `p_d'
			gen num=$n
			keep parm p num
			duplicates drop
			save "$gsdRawTemp\1qvalue_ate_difference3_psy_$n", replace
			restore
	local replace
	local append append
	global n=$n+1
}

		foreach k of numlist 1/4 {
		  matrix stars`k'[1,1] = (abs(A`k'[1,2]) < 0.1) + (abs(A`k'[1,2]) < 0.05) + (abs(A`k'[1,2]) < 0.01)
		  matrix stars`k'[1,4] = (abs(A`k'[1,5]) < 0.1) + (abs(A`k'[1,5]) < 0.05) + (abs(A`k'[1,5]) < 0.01)
		  matrix stars`k'[1,7] = (abs(A`k'[1,8]) < 0.1) + (abs(A`k'[1,8]) < 0.05) + (abs(A`k'[1,8]) < 0.01)
		  matrix stars`k'[1,10] = (abs(A`k'[1,11]) < 0.1) + (abs(A`k'[1,11]) < 0.05) + (abs(A`k'[1,11]) < 0.01)
		  matrix stars`k'[1,13] = (abs(A`k'[1,14]) < 0.1) + (abs(A`k'[1,14]) < 0.05) + (abs(A`k'[1,14]) < 0.01)
		  matrix stars`k'[1,16] = (abs(A`k'[1,17]) < 0.1) + (abs(A`k'[1,17]) < 0.05) + (abs(A`k'[1,17]) < 0.01)
		  }
		  
		foreach k of numlist 1/5 {
		  matrix starZ`k'[1,1] = (abs(O`k'[1,2]) < 0.1) + (abs(O`k'[1,2]) < 0.05) + (abs(O`k'[1,2]) < 0.01)
		  matrix starZ`k'[1,4] = (abs(O`k'[1,5]) < 0.1) + (abs(O`k'[1,5]) < 0.05) + (abs(O`k'[1,5]) < 0.01)
		  matrix starZ`k'[1,7] = (abs(O`k'[1,8]) < 0.1) + (abs(O`k'[1,8]) < 0.05) + (abs(O`k'[1,8]) < 0.01)
		  matrix starZ`k'[1,10] = (abs(O`k'[1,11]) < 0.1) + (abs(O`k'[1,11]) < 0.05) + (abs(O`k'[1,11]) < 0.01)
		  matrix starZ`k'[1,13] = (abs(O`k'[1,14]) < 0.1) + (abs(O`k'[1,14]) < 0.05) + (abs(O`k'[1,14]) < 0.01)
		  matrix starZ`k'[1,16] = (abs(O`k'[1,17]) < 0.1) + (abs(O`k'[1,17]) < 0.05) + (abs(O`k'[1,17]) < 0.01)
		  }
		
		
	
preserve
do "${gsdDo}\3-2-a-AdjustedPValues.do"
restore
	
	foreach num of numlist 1/4{
		mat A`num'[1,3] = B3[`num',1]
		}
	foreach num of numlist 1/4{
		mat A`num'[1,6] = B4[`num',1]
		}
	foreach num of numlist 1/4{
		mat A`num'[1,9] = B8[`num',1]
		}
	foreach num of numlist 1/4{
		mat A`num'[1,12] = B5[`num',1]
		}
	foreach num of numlist 1/4{
		mat A`num'[1,15] = B6[`num',1]
		}
	foreach num of numlist 1/4{
		mat A`num'[1,18] = B9[`num',1]
		}

preserve
do "${gsdDo}\3-2-b-AdjustedPValuesPsy.do"
restore
	foreach num of numlist 1/5{
		mat O`num'[1,3] = C3[`num',1]
		}
	foreach num of numlist 1/5{
		mat O`num'[1,6] = C4[`num',1]
		}
	foreach num of numlist 1/5{
		mat O`num'[1,9] = C8[`num',1]
		}
	foreach num of numlist 1/5{
		mat O`num'[1,12] = C5[`num',1]
		}
	foreach num of numlist 1/5{
		mat O`num'[1,15] = C6[`num',1]
		}
	foreach num of numlist 1/5{
		mat O`num'[1,18] = C9[`num',1]
		}

********************************************************************************

* Make table

		matrix rownames A1 = "Employment index" 
		frmttable , statmat(A1) replace sdec(3)  substat(2) plain annotate(stars1) asymbol(*,**,***)   /// 
		title (Table 4: Local average treatment effects on main outcomes) ///
		ctitles("") 
		frmttable , statmat(A2) replace sdec(3) substat(2) annotate(stars2) asymbol(*,**,***)  append rtitles("Consumption index") 
		frmttable , statmat(A3) replace sdec(3) substat(2) annotate(stars3) asymbol(*,**,***)  append rtitles("Savings index") 
		frmttable , statmat(A4) replace sdec(3) substat(2) annotate(stars4) asymbol(*,**,***)  append rtitles("Business skills index") 
		frmttable , statmat(O1) replace sdec(3) substat(2) annotate(starZ1) asymbol(*,**,***)  append rtitles("Psychological index") 
		frmttable , statmat(O2) replace sdec(3) substat(2) annotate(starZ2) asymbol(*,**,***)  append rtitles("Risk index") 
		frmttable , statmat(O3) replace sdec(3) substat(2) annotate(starZ3) asymbol(*,**,***)  append rtitles("Trust index") 
		frmttable , statmat(O4) replace sdec(3) substat(2) annotate(starZ4) asymbol(*,**,***)  append rtitles("Crime index") 
		frmttable , statmat(O5) replace sdec(3) substat(2) annotate(starZ5) asymbol(*,**,***)  append rtitles("Migration index") 
		frmttable , statmat(N) replace sdec(0) append rtitles("Observations") 
		frmttable using "$gsdTables/Table4", statmat(D) replay replace tex sdec(3) append rtitles("F-stat") plain
		